Take-home Exercise 2

Spatio-temporal Analysis of COVID-19 Vaccination Trends at the Sub-district Level, DKI Jakarta

Published

February 26, 2023

Modified

March 1, 2023

1 Context

1.1 Problem Statement

Understanding how Covid-19 vaccination rate is distributed around DKI Jakarta and how it changes over time.

1.2 Data

Data Format (Type) Description Source
Riwayat File Vaksinasi DKI Jakarta Excel (Aspatial) Details the daily number of vaccinations in Jakarta by person type and sub-district (Only the first day of each month is used) Riwatat File Vaksinasi DKI Jakarta (2022)
DKI Jakarta Adminstration Boundary 2019 Shapefile (Geospatial) DKI Jakarta administrative boundaries Indonesia Geospasial (2019)
Note

For the Vaksinasi DKI Jakarta data for March 2022, there was no available link to 1 March 2022 data. Thus for March 2022, the first day available, 2 March 2022, will be used instead.

2 Setting Up

2.1 Installing and Loading Packages

pacman::p_load(sf, tmap, sfdep, tidyverse, zoo, readxl, plotly, Kendall, ggplot2, dplyr)

  • sf: importing and handling geospatial data

  • tidyverse: wrangling attribute data

  • sfdep: newer package of spdep (does not require conversion to sp)

  • tmap: prepare cartographic quality chropleth maps

  • zoo

  • readxl: read excel (xlsx) format

  • plotly:

  • Kendall: to perform Mann-Kendall test

  • ggplot2: plot temporal trend

  • dplyr: join_by() function

2.2 DKI Jakarta Administrative Boundary 2019

2.2.1 Importing Data

To import a shapefile, use st_read() where dsn is the targeted directory and layer is the name of the files.

jkt <- st_read(
  dsn = "data/geospatial",
  layer = "BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA"
)
Reading layer `BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA' from data source 
  `D:\bellekwang\IS415\take_home_ex\THE2\data\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 269 features and 161 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 106.3831 ymin: -6.370815 xmax: 106.9728 ymax: -5.184322
Geodetic CRS:  WGS 84
  • jkt is a multipolygon sf dataframe with 269 features

2.2.2 Projection

The current geodetic CRS is WGS 84 which is wrong. Using st_transform(), change the CRS from WGS 84 to DGN95 (EPSG:23845), the national projected coordinates systems of Indonesia TM-3 zone 54.1.

jkt <- st_transform(jkt, crs = 23845)
st_crs(jkt)
Coordinate Reference System:
  User input: EPSG:23845 
  wkt:
PROJCRS["DGN95 / Indonesia TM-3 zone 54.1",
    BASEGEOGCRS["DGN95",
        DATUM["Datum Geodesi Nasional 1995",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4755]],
    CONVERSION["Indonesia TM-3 zone 54.1",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",139.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9999,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",200000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",1500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre."],
        AREA["Indonesia - onshore east of 138°E."],
        BBOX[-9.19,138,-1.49,141.01]],
    ID["EPSG",23845]]

2.2.3 Data Cleaning

2.2.3.1 Outer Island

First, using tmap functions tm_shape() and tm_polygons(), plot out jkt to understand the area covered. I will also be adding the “DESA” label to identify which are the areas to exclude.

tmap_mode("view")
tm_shape(jkt) +
  tm_polygons() +
  tm_text("DESA") +
tm_view(set.zoom.limits = c(9,12))

There are many outer islands shown on the plot above. As seen from the plot, the areas to be excluded are: Pulau Harapan, Pulau Kelapa, Pulau Panggang, Pulau Tidung, Pulau Pari and Pulau Untung Jawa.

outer_islands <- c("PULAU HARAPAN", "PULAU KELAPA", "PULAU PANGGANG", "PULAU TIDUNG", "PULAU PARI", "PULAU UNTUNG JAWA")

Use the filter() function to remove the outer islands.

jkt <- jkt %>%
  filter(!DESA %in% outer_islands)

Now using the plot mode, plot out jkt and ensure all outer islands are excluded.

tmap_mode("plot")
tm_shape(jkt) +
  tm_polygons()

As you can see from the plot, all the outer islands are excluded.

2.2.3.2 Unnecessary Fields

Only the first nine columns of the sf dataframe is relevant. Using the select() function, pick out the first nine columns.

jkt <- select(jkt, 1:9)

2.2.4 Final jkt Data Frame

st_geometry(jkt)
Geometry set for 263 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -3644275 ymin: 663887.8 xmax: -3606237 ymax: 701380.1
Projected CRS: DGN95 / Indonesia TM-3 zone 54.1
First 5 geometries:

2.3 Riwayat File Vaksinasi DKI Jakarta

There are 12 excel files from July 2021 to June 2022 to be imported. Using read_excel() function, import the “Data Kelurahan” sheet.

july_2021 = read_excel("data/aspatial/july_2021.xlsx", sheet = 1)
aug_2021 = read_excel("data/aspatial/aug_2021.xlsx", sheet = "Data Kelurahan")
sept_2021 = read_excel("data/aspatial/sept_2021.xlsx", sheet = "Data Kelurahan")
oct_2021 = read_excel("data/aspatial/oct_2021.xlsx", sheet = "Data Kelurahan")
nov_2021 = read_excel("data/aspatial/nov_2021.xlsx", sheet = "Data Kelurahan")
dec_2021 = read_excel("data/aspatial/dec_2021.xlsx", sheet = "Data Kelurahan")
jan_2022 = read_excel("data/aspatial/jan_2022.xlsx", sheet = "Data Kelurahan")
feb_2022 = read_excel("data/aspatial/feb_2022.xlsx", sheet = "Data Kelurahan")
mar_2022 = read_excel("data/aspatial/mar_2022.xlsx", sheet = "Data Kelurahan")
apr_2022 = read_excel("data/aspatial/apr_2022.xlsx", sheet = "Data Kelurahan")
may_2022 = read_excel("data/aspatial/may_2022.xlsx", sheet = "Data Kelurahan")
june_2022 = read_excel("data/aspatial/june_2022.xlsx", sheet = "Data Kelurahan")

2.3.1 Data Cleaning

I will exclude all the outer_islands identified earlier.

july_2021 <- july_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
aug_2021 <- aug_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
sept_2021 <- sept_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
oct_2021 <- oct_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
nov_2021 <- nov_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
dec_2021 <- dec_2021 %>%
  filter(!KELURAHAN %in% outer_islands)
jan_2022 <- jan_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
feb_2022 <- feb_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
mar_2022 <- mar_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
apr_2022 <- apr_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
may_2022 <- may_2022 %>%
  filter(!KELURAHAN %in% outer_islands)
june_2022 <- june_2022 %>%
  filter(!KELURAHAN %in% outer_islands)

All the tables now have 262 observations, 1 less than the jkt multipolygon dataframe.

2.3.1.1 Unnecessary Columns

Only the first 6 columns are important to the exercise.

july_2021 <- select(july_2021, 1:6)
aug_2021 <- select(aug_2021, 1:6)
sept_2021 <- select(sept_2021, 1:6)
oct_2021 <- select(oct_2021, 1:6)
nov_2021 <- select(nov_2021, 1:6)
dec_2021 <- select(dec_2021, 1:6)
jan_2022 <- select(jan_2022, 1:6)
feb_2022 <- select(feb_2022, 1:6)
mar_2022 <- select(mar_2022, 1:6)
apr_2022 <- select(apr_2022, 1:6)
may_2022 <- select(may_2022, 1:6)
june_2022 <- select(june_2022, 1:6)

Now all the datasets only have 6 columns.

2.3.2 Combining Datasets

First, add a date column to each data frame using the mutate() and as.Date() function.

july_2021 <- july_2021 %>%
  mutate(dateMonth = as.Date("2021-07-01"))
aug_2021 <- aug_2021 %>%
  mutate(dateMonth = as.Date("2021-08-01"))
sept_2021 <- sept_2021 %>%
  mutate(dateMonth = as.Date("2021-09-01"))
oct_2021 <- oct_2021 %>%
  mutate(dateMonth = as.Date("2021-10-01"))
nov_2021 <- nov_2021 %>%
  mutate(dateMonth = as.Date("2021-11-01"))
dec_2021 <- dec_2021 %>%
  mutate(dateMonth = as.Date("2021-12-01"))
jan_2022 <- jan_2022 %>%
  mutate(dateMonth = as.Date("2022-01-01"))
feb_2022 <- feb_2022 %>%
  mutate(dateMonth = as.Date("2022-02-01"))
mar_2022 <- mar_2022 %>%
  mutate(dateMonth = as.Date("2022-03-02"))
apr_2022 <- apr_2022 %>%
  mutate(dateMonth = as.Date("2022-04-01"))
may_2022 <- may_2022 %>%
  mutate(dateMonth = as.Date("2022-05-01"))
june_2022 <- june_2022 %>%
  mutate(dateMonth = as.Date("2022-06-01"))
july_2021 <- july_2021 %>%
  mutate(dateMonthNum = as.numeric(20210701))
aug_2021 <- aug_2021 %>%
  mutate(dateMonthNum = as.numeric(20210801))
sept_2021 <- sept_2021 %>%
  mutate(dateMonthNum = as.numeric(20210901))
oct_2021 <- oct_2021 %>%
  mutate(dateMonthNum = as.numeric(20211001))
nov_2021 <- nov_2021 %>%
  mutate(dateMonthNum = as.numeric(20211101))
dec_2021 <- dec_2021 %>%
  mutate(dateMonthNum = as.numeric(20211201))
jan_2022 <- jan_2022 %>%
  mutate(dateMonthNum = as.numeric(20220101))
feb_2022 <- feb_2022 %>%
  mutate(dateMonthNum = as.numeric(20220201))
mar_2022 <- mar_2022 %>%
  mutate(dateMonthNum = as.numeric(20220301))
apr_2022 <- apr_2022 %>%
  mutate(dateMonthNum = as.numeric(20220401))
may_2022 <- may_2022 %>%
  mutate(dateMonthNum = as.numeric(20220501))
june_2022 <- june_2022 %>%
  mutate(dateMonthNum = as.numeric(20220601))

Next, I will combine all the 12 datasets into 1 dataset using rbind() function.

vacc_data <- rbind(july_2021, aug_2021, sept_2021, oct_2021, nov_2021, dec_2021, jan_2022, feb_2022, mar_2022, apr_2022, may_2022, june_2022)

2.3.3 Final Dataset

head(vacc_data)
# A tibble: 6 × 8
  `KODE KELURAHAN` WILAYAH …¹ KECAM…² KELUR…³ SASARAN BELUM…⁴ dateMonth  dateM…⁵
  <chr>            <chr>      <chr>   <chr>     <dbl>   <dbl> <date>       <dbl>
1 <NA>             <NA>       <NA>    TOTAL   7739060 5041111 2021-07-01  2.02e7
2 3172051003       JAKARTA U… PADEMA… ANCOL     20393   13272 2021-07-01  2.02e7
3 3173041007       JAKARTA B… TAMBORA ANGKE     25785   16477 2021-07-01  2.02e7
4 3175041005       JAKARTA T… KRAMAT… BALE K…   25158   18849 2021-07-01  2.02e7
5 3175031003       JAKARTA T… JATINE… BALI M…    8683    5743 2021-07-01  2.02e7
6 3175101006       JAKARTA T… CIPAYU… BAMBU …   22768   15407 2021-07-01  2.02e7
# … with abbreviated variable names ¹​`WILAYAH KOTA`, ²​KECAMATAN, ³​KELURAHAN,
#   ⁴​`BELUM VAKSIN`, ⁵​dateMonthNum

The final vacc_data dataset has a total of 7 columns:

  1. Identifiers: Kode Kelurahan / Kelurahan
  2. Targeted people to be vaccinated: Sasaran
  3. Yet to be vaccinated: Belum Vaksin
  4. Date identifier: dateMonth

3 Choropleth Mapping and Analysis

3.1 Monthly Vaccination Rates at Sub-district Level

3.1.1 Computing Monthly Vaccination Rates

  • SASARAN (Target - target people to be vaccinated)

  • BELUM VAKSIN (Yet to be vaccinated)

Using the equation below, compute the monthly vaccination rate and add the rate value as a new column ‘monthly_vacc_rate’.

\[ Monthly Vaccination Rate = \frac{(SASARAN - BELUMVAKSIN)}{SUM(SASARAN, BELUMVAKSIN)} \]

vacc_data <- vacc_data %>%
  mutate(monthly_vacc_rate = (SASARAN - `BELUM VAKSIN`)/(SASARAN + `BELUM VAKSIN`))

3.1.2 Combining Dataset and Shapefile

Using the left_join() function, with identifiers “DESA” and “KELURAHAN”, combine the vacc_data dataset with the jkt shapefile.

vacc_july_2021 <- vacc_data %>%
  filter(dateMonth == "2021-07-01")
jkt_july_2021 <- left_join(jkt, vacc_july_2021, by = c("DESA" = "KELURAHAN" ))

3.2 Choropleth Mapping

Using the qtm() function, plot out the choropleth map for July 2021.

qtm(jkt_july_2021, "monthly_vacc_rate", title = "July 2021", fill.palette="Greens")

As seen above, there are missing monthly vaccinated rate values for some sub-districts. There is a mismatch and lack of vaccination data for these areas.

I will leave them as missing values instead of giving a monthly vaccination rate of 0 to avoid confusion with actual areas with recorded 0 rates.

3.2.1 Other Months

Repeat the above steps for the next 11 months.

Code 
#AUGUST 2021
vacc_aug_2021 <- vacc_data %>%
  filter(dateMonth == "2021-08-01")
jkt_aug_2021 <- left_join(jkt, vacc_aug_2021, by = c("DESA" = "KELURAHAN" ))

#SEPTEMBER 2021
vacc_sept_2021 <- vacc_data %>%
  filter(dateMonth == "2021-09-01")
jkt_sept_2021 <- left_join(jkt, vacc_sept_2021, by = c("DESA" = "KELURAHAN" ))

#OCTOBER 2021
vacc_oct_2021 <- vacc_data %>%
  filter(dateMonth == "2021-10-01")
jkt_oct_2021 <- left_join(jkt, vacc_oct_2021, by = c("DESA" = "KELURAHAN" ))

#NOVEMBER 2021
vacc_nov_2021 <- vacc_data %>%
  filter(dateMonth == "2021-11-01")
jkt_nov_2021 <- left_join(jkt, vacc_nov_2021, by = c("DESA" = "KELURAHAN" ))

#DECEMBER 2021
vacc_dec_2021 <- vacc_data %>%
  filter(dateMonth == "2021-12-01")
jkt_dec_2021 <- left_join(jkt, vacc_dec_2021, by = c("DESA" = "KELURAHAN" ))

#JANUARY 2022
vacc_jan_2022 <- vacc_data %>%
  filter(dateMonth == "2022-01-01")
jkt_jan_2022 <- left_join(jkt, vacc_jan_2022, by = c("DESA" = "KELURAHAN" ))

#FEBRUARY 2022
vacc_feb_2022 <- vacc_data %>%
  filter(dateMonth == "2022-02-01")
jkt_feb_2022 <- left_join(jkt, vacc_feb_2022, by = c("DESA" = "KELURAHAN" ))

#MARCH 2022
vacc_mar_2022 <- vacc_data %>%
  filter(dateMonth == "2022-03-02")
jkt_mar_2022 <- left_join(jkt, vacc_mar_2022, by = c("DESA" = "KELURAHAN" ))

#APRIL 2022
vacc_apr_2022 <- vacc_data %>%
  filter(dateMonth == "2022-04-01")
jkt_apr_2022 <- left_join(jkt, vacc_apr_2022, by = c("DESA" = "KELURAHAN" ))

#MAY 2022
vacc_may_2022 <- vacc_data %>%
  filter(dateMonth == "2022-05-01")
jkt_may_2022 <- left_join(jkt, vacc_may_2022, by = c("DESA" = "KELURAHAN" ))

#JUNE 2022
vacc_june_2022 <- vacc_data %>%
  filter(dateMonth == "2022-06-01")
jkt_june_2022 <- left_join(jkt, vacc_june_2022, by = c("DESA" = "KELURAHAN" ))
Code 
qtm(jkt_july_2021, "monthly_vacc_rate", title = "July 2021", fill.palette="Greens")

Code 
qtm(jkt_aug_2021, "monthly_vacc_rate", title = "Aug 2021", fill.palette="Greens")

Code 
qtm(jkt_sept_2021, "monthly_vacc_rate", title = "Sept 2021", fill.palette="Greens")

Code 
qtm(jkt_oct_2021, "monthly_vacc_rate", title = "Oct 2021", fill.palette="Greens")

Code 
qtm(jkt_nov_2021, "monthly_vacc_rate", title = "Nov 2021", fill.palette="Greens")

Code 
qtm(jkt_dec_2021, "monthly_vacc_rate", title = "Dec 2021", fill.palette="Greens")

Code 
qtm(jkt_jan_2022, "monthly_vacc_rate", title = "Jan 2022", fill.palette="Greens")

Code 
qtm(jkt_feb_2022, "monthly_vacc_rate", title = "Feb 2022", fill.palette="Greens")

Code 
qtm(jkt_mar_2022, "monthly_vacc_rate", title = "Mar 2022", fill.palette="Greens")

Code 
qtm(jkt_apr_2022, "monthly_vacc_rate", title = "Apr 2022", fill.palette="Greens")

Code 
qtm(jkt_may_2022, "monthly_vacc_rate", title = "May 2022", fill.palette="Greens")

Code 
qtm(jkt_june_2022, "monthly_vacc_rate", title = "June 2022", fill.palette="Greens")

3.3 Spatial Patterns

4 Local Gi* Analysis

4.1 Computing Local Gi* of Monthly Vaccination Rate

To compute local Gi* statistics, it requires the neighbour list (using st_contiguity()) and the wts (using st_inverse_distance()).

Using jkt_july_2021 as an example:

jkt_july_2021 = jkt_july_2021 %>% drop_na(monthly_vacc_rate)
jkt_july_2021_idw <- jkt_july_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

Next, add in the local Gi* values into a new sf dataset HCSA_july_2021 as a new column with the mutate() function, adding .before = 1 to place the new column as the first column. The unnest() function is used to this make each element of the local_Gi list its own row

HCSA_july_2021 <- jkt_july_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)
HCSA_july_2021
Simple feature collection with 252 features and 27 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -3644275 ymin: 663887.8 xmax: -3606237 ymax: 701380.1
Projected CRS: DGN95 / Indonesia TM-3 zone 54.1
# A tibble: 252 × 28
   gi_star    e_gi      var_gi p_value p_sim p_fol…¹ skewn…² kurto…³ nb    wts  
     <dbl>   <dbl>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>   <dbl> <nb>  <lis>
 1  1.81   0.00396     1.26e-7 0.0707   0.08    0.04  0.291  -0.360  <int> <dbl>
 2  2.79   0.00399     9.26e-8 0.00524  0.04    0.02  0.669   0.179  <int> <dbl>
 3 -0.544  0.00401     1.36e-7 0.586    0.52    0.26  0.373   0.165  <int> <dbl>
 4 -0.928  0.00395     9.61e-8 0.353    0.3     0.15  0.460   0.220  <int> <dbl>
 5  2.35   0.00394     6.89e-8 0.0186   0.04    0.02 -0.0269  0.0678 <int> <dbl>
 6  1.21   0.00397     9.30e-8 0.226    0.24    0.12  0.289  -0.0724 <int> <dbl>
 7  0.970  0.00399     5.41e-8 0.332    0.28    0.14  0.395   0.313  <int> <dbl>
 8  0.0232 0.00395     1.48e-7 0.982    0.96    0.48  0.172  -0.121  <int> <dbl>
 9  0.821  0.00400     5.34e-8 0.411    0.46    0.23  0.110  -0.623  <int> <dbl>
10 -0.330  0.00398     7.74e-8 0.742    0.7     0.35  0.148  -0.317  <int> <dbl>
# … with 242 more rows, 18 more variables: OBJECT_ID <dbl>, KODE_DESA <chr>,
#   DESA <chr>, KODE <dbl>, PROVINSI <chr>, KAB_KOTA <chr>, KECAMATAN.x <chr>,
#   DESA_KELUR <chr>, JUMLAH_PEN <dbl>, `KODE KELURAHAN` <chr>,
#   `WILAYAH KOTA` <chr>, KECAMATAN.y <chr>, SASARAN <dbl>,
#   `BELUM VAKSIN` <dbl>, dateMonth <date>, dateMonthNum <dbl>,
#   monthly_vacc_rate <dbl>, geometry <MULTIPOLYGON [m]>, and abbreviated
#   variable names ¹​p_folded_sim, ²​skewness, ³​kurtosis

4.2 Displaying Gi* Maps

Using the plot mode, plot out both the Gi* Map and the p-value of Gi* using the HCDA_july_2021 sf dataset containing the newly added Gi* values.

tmap_mode("plot")
map1 <- tm_shape(HCSA_july_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of July 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_july_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3 Other Months

Repeat the previous steps to the other 11 months.

4.3.1 Aug 2021

Code 
jkt_aug_2021 = jkt_aug_2021 %>% drop_na(monthly_vacc_rate)
jkt_aug_2021_idw <- jkt_aug_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_aug_2021 <- jkt_aug_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_aug_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of August 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_aug_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.2 September 2021

Code 
jkt_sept_2021 = jkt_sept_2021 %>% drop_na(monthly_vacc_rate)
jkt_sept_2021_idw <- jkt_sept_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_sept_2021 <- jkt_sept_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_sept_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of September 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_sept_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.3 October 2021

Code 
jkt_oct_2021 = jkt_oct_2021 %>% drop_na(monthly_vacc_rate)
jkt_oct_2021_idw <- jkt_oct_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_oct_2021 <- jkt_oct_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_oct_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of October 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_oct_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.4 November 2022

Code 
jkt_nov_2021 = jkt_nov_2021 %>% drop_na(monthly_vacc_rate)
jkt_nov_2021_idw <- jkt_nov_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_nov_2021 <- jkt_nov_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_nov_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of November 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_nov_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.5 December 2021

Code 
jkt_dec_2021 = jkt_dec_2021 %>% drop_na(monthly_vacc_rate)
jkt_dec_2021_idw <- jkt_dec_2021 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_dec_2021 <- jkt_dec_2021_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_dec_2021) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of December 2021",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_dec_2021) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.6 January 2022

Code 
jkt_jan_2022 = jkt_jan_2022 %>% drop_na(monthly_vacc_rate)
jkt_jan_2022_idw <- jkt_jan_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_jan_2022 <- jkt_jan_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_jan_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of January 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_jan_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.7 February 2022

Code 
jkt_feb_2022 = jkt_feb_2022 %>% drop_na(monthly_vacc_rate)
jkt_feb_2022_idw <- jkt_feb_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_feb_2022 <- jkt_feb_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_feb_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of February 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_feb_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.8 March 2022

Code 
jkt_mar_2022 = jkt_mar_2022 %>% drop_na(monthly_vacc_rate)
jkt_mar_2022_idw <- jkt_mar_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_mar_2022 <- jkt_mar_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_mar_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of March 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_mar_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.9 April 2022

Code 
jkt_apr_2022 = jkt_apr_2022 %>% drop_na(monthly_vacc_rate)
jkt_apr_2022_idw <- jkt_apr_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_apr_2022 <- jkt_apr_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_apr_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of April 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_apr_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.10 May 2022

Code 
jkt_may_2022 = jkt_may_2022 %>% drop_na(monthly_vacc_rate)
jkt_may_2022_idw <- jkt_may_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_may_2022 <- jkt_may_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_may_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of May 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_may_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.3.11 June 2022

Code 
jkt_june_2022 = jkt_june_2022 %>% drop_na(monthly_vacc_rate)
jkt_june_2022_idw <- jkt_june_2022 %>%
  mutate(
    nb = st_contiguity(geometry),
    wts = st_inverse_distance(
      nb, geometry, scale = 1, alpha = 1
    ),
  .before = 1
  )

HCSA_june_2022 <- jkt_june_2022_idw %>% 
  mutate(local_Gi = local_gstar_perm(
    monthly_vacc_rate, nb, wt, nsim = 99),
         .before = 1) %>%
  unnest(local_Gi)

tmap_mode("plot")
map1 <- tm_shape(HCSA_june_2022) +
  tm_fill("gi_star", palette = "-RdBu") + 
  tm_borders(alpha = 0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Gi* of June 2022",
            main.title.size = 0.8)

map2 <- tm_shape(HCSA_june_2022) +
  tm_fill("p_value",
          breaks = c(0, 0.001, 0.01, 0.05, 1),
              labels = c("0.001", "0.01", "0.05", "Not sig"), palette = "-Greens") + 
  tm_borders(alpha = 0.5) +
  tm_layout(main.title = "p-value of Gi*",
            main.title.size = 0.8)

tmap_arrange(map1, map2, ncol = 2)

4.4 Statistical Conclusion

5 Emerging Hot Spot Analysis

5.1 Mann-Kendall Test

(using spatio-temporal local Gi* values)

5.1.1 Creating Time Series Cube

colnames(jkt)[3] <- "KELURAHAN"
#vacc_data = vacc_data %>% drop_na(c("KODE KELURAHAN", "WILAYAH KOTA", "KECAMATAN"))
#vacc_data_st <- spacetime(vacc_data, jkt, ".KELURAHAN", "dateMonth")
vacc_data_st <- spacetime(
  .data = vacc_data, 
  .geometry = jkt, 
  .loc_col = "KELURAHAN", 
  .time_col = "dateMonthNum")
vacc_data_st_complete <- complete_spacetime_cube(vacc_data_st)
is_spacetime_cube(vacc_data_st_complete)
[1] TRUE

5.1.2 Computing Spatial Weights

vacc_st_nb <- vacc_data_st_complete %>%
  activate("geometry") %>%
  mutate(nb = include_self(st_contiguity(geometry)),
         wt = st_inverse_distance(nb, geometry,
                                  scale = 1,
                                  alpha = 1),
         .before = 1) %>%
  set_nbs("nb") %>%
  set_wts("wt")
head(vacc_st_nb)
# A tibble: 6 × 11
  KELURAHAN   dateM…¹ KODE …² WILAY…³ KECAM…⁴ SASARAN BELUM…⁵ dateMonth  month…⁶
  <chr>         <dbl> <chr>   <chr>   <chr>     <dbl>   <dbl> <date>       <dbl>
1 KEAGUNGAN    2.02e7 317303… JAKART… TAMAN …   15262   10275 2021-07-01   0.195
2 GLODOK       2.02e7 317303… JAKART… TAMAN …    7192    3561 2021-07-01   0.338
3 HARAPAN MU…  2.02e7 317103… JAKART… KEMAYO…   19987   13623 2021-07-01   0.189
4 CEMPAKA BA…  2.02e7 317103… JAKART… KEMAYO…   27330   18417 2021-07-01   0.195
5 PASAR BARU   2.02e7 317102… JAKART… SAWAH …   11656    6067 2021-07-01   0.315
6 KARANG ANY…  2.02e7 317102… JAKART… SAWAH …   23056   15040 2021-07-01   0.210
# … with 2 more variables: nb <list>, wt <list>, and abbreviated variable names
#   ¹​dateMonthNum, ²​`KODE KELURAHAN`, ³​`WILAYAH KOTA`, ⁴​KECAMATAN,
#   ⁵​`BELUM VAKSIN`, ⁶​monthly_vacc_rate

5.1.3 Computing Gi*

length(vacc_st_nb$wt)

vacc_st_nb <- vacc_st_nb %>% drop_na("monthly_vacc_rate")
gi_stars <- vacc_st_nb %>% 
  mutate(gi_star = local_gstar_perm(
    monthly_vacc_rate, nb, wt)) %>% 
  unnest(gi_star)

5.1.4 Mann-Kendall Test

ehsa <- gi_stars %>%
  group_by(KELURAHAN) %>%
  summarise(mk = list(
    unclass(MannKendall(gi_star)))) %>%
  tidyr::unnest_wider(mk)

5.3 EHSA Map

(only significant values)

vacc_data_st_complete$monthly_vacc_rate[is.na(vacc_data_st_complete$monthly_vacc_rate)] <- 0
ehsa <- emerging_hotspot_analysis(
  x = vacc_data_st_complete, 
  .var = "monthly_vacc_rate", 
  k = 1, 
  nsim = 99
)
jkt_ehsa <- jkt %>%
  left_join(ehsa, by = c("KELURAHAN" = "location" ))
ehsa_sig <- jkt_ehsa  %>%
  filter(p_value < 0.05)
tmap_mode("plot")
tm_shape(jkt_ehsa) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
tm_shape(ehsa_sig) +
  tm_fill("classification") + 
  tm_borders(alpha = 0.4)

5.4 Spatial Patterns

6 Conclusion